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Abstract 

The Minimum Ignition Energy (MIE) of an initially Gaussian temperature 
profile is found both by Direct Numerical Simulations (DNS) and from a new 
novel model. The model is based on solving the heat diffusion equation in 
zero dimensions for a Gaussian velocity distribution. The chemistry is taken 
into account through the ignition delay time, which is required as input to 
the model. The model results reproduce the DNS results very well for the 
Hydrogen mixture investigated. 

Furthermore, the effect of ignition source dimensionality is explored, and 
it is shown that for compact ignition kernels there is a strong effect on dimen- 
sionality. Here, three, two and one dimensional ignition sources represent a 
spherical kernel, a long spark and an ignition sheet, respectively. 
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1. Introduction 

A flammable mixture of fuel and oxidant exposed to a flame will react 
and produce heat and combustion products in a very short time. If the same 
mixture is not exposed to a significant heat source at any time, the fuel- 
oxidant mixture is, however, stable and will not react. For e.g. safety issues 
it is important to know the amount of energy required in order to ignite the 
mixture and initiate the combustion process. It turns out that in addition to 
temperature and pressure, the Minimum Ignition Energy (MIE) for a given 
mixture depends on at least three different parameters; the geometry of the 



Preprint submitted to Combustion and Flame 



October 7, 2011 



ignition source, the radius of the ignition source r s and the deposition period 
td- For deposition periods in the range t a « td « t c , the MIE should 
be independent on td when t a = Sf/c is the acoustic time scale, t c = ajS\ 
is the chemical time scale, 5f ~ ol/Sl is the thickness of the flame, a is 
the thermal diffusivity, Sl is the laminar flame speed and c is the speed of 
sound. Several author groups calculated the MIE using different numerical 
techniques [H EJ [3j H] and experimental investigations [5]. 

Depending on the dimensionality of the heat source the MIE will vary 
considerably. The initial ignition kernel may have any given profile, depend- 
ing on the heat source. In the current work a Gaussian profile has been 
chosen. A one dimensional heat source correspond to heating the mixture 
in an infinitely large plane, which could possibly be realized by laser sheets. 
A two dimensional heat source correspond to an infinitely long cylindrical 
ignition kernel. This is in essence the same as a spark ignition where the 
length of the spark is significantly longer than its width. Finally a Gaussian 
profile in three dimensions correspond to a spherical source which could be 
thought of as a very short spark, or as ignition by the combined effect of 
several focused lasers. 

In the present work it is assumed that all thermal energy is deposited 
with a Gaussian distribution and constant pressure prior to starting the cal- 
culation. This essentially means that t a << td << t c . This simplification 
has been chosen in order to remove one variable from the equation and conse- 
quently to more easily see the fundamental underlying physics. Furthermore 
the focus is on Hydrogen, but the methods and conclusions should qualita- 
tively be independent on the fuel and therefore be of generic interest. 

2. Numerical simulations 

The minimum ignition energy for a combustible hydrogen mixture is 
found be running a series of Direct Numerical Simulations (DNS) with the 
Pencil Code[6l [7j. Unlike Large Eddy Simulations (LES) and Reynold- 
average Navier-Stokes simulations (RANS), which use turbulence modelling, 
DNS solves the full Navier-Stokes equations without the use of any modelling 
and filtering. 

To simulate the physical problem the conservation equations for mass, 
momentum, species and energy have to be solved together with the equation 
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of state. The equation for conservation of mass is given as 

^ = -v-u, (1) 

where U is the velocity vector, p is the density and D/Dt = d/dt + U ■ V is 
the advective derivative. 

The momentum equation has the form 

^ = i("Vp + F„), (2) 

where F„ s = V ■ (2pvS) is the viscous force, S is the rate of strain tensor, 
and p is the pressure. The species evolution equation is 

^ = ^(-V-J,+u),), (3) 

where uj is the reaction rate, Jj is the diffusive flux and Yj is the mass fraction 
of specie j. Lastly, the energy equation is solved for the temperature 

R\DlnT =sr DYj(R hA R 2vS* V-q 

' M J Dt j Dt \ M i T J M T pT ' 1 ] 

where c p is the heat capacity at constant pressure, R is the universal gas 
constant, M is the molar mass, T is the temperature, h is the enthalpy and 
q is the heat flux. For temporal discretization the Pencil Code uses a third 
order Runge-Kutta method, while a sixth-order central difference scheme is 
used for the spatial discretization. For a more thorough discussion on the 
equations solved see Babkovskaia, Haugen and Brandenburg (2011) [7]. 

To simulate an ignition source an initial temperature distribution is im- 
posed in the domain. The distribution is chosen to be Gaussian, and parallels 
could be drawn to a real life heat source that has its hottest spot in the center. 
For instance, if a cross-section of a real electrical spark is made, the temper- 
ature distribution at this cross-section could compare well with a Gaussian 
shaped temperature distribution in two dimensions. The initial distribution 
is given by 

T(r) = (T max -T )e"fe) + T , (5) 

where T max is the maximum temperature, T is ambient temperature and r so 
is the radius of the distribution. Fig. [T] illustrates the initial temperature 
profile in one dimension for one of the simulations. 
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Figure 1: Illustration of the initial Gaussian temperature distribution, where T = 300 K, 
r SQ = 0.05 cm, T max = 1300 K and the size of the domain L — 0.3 cm. 

To simulate a closed vessel or container the spatial derivative of the tem- 
perature, species and the density, together with the value of all velocity 
components, are set to zero at the walls. The domain is chosen to be suffi- 
ciently large, which implies that all spatial gradients are close to zero at the 
boundaries for all relevant times. The size of the domain, L, ranges from 
0.5 cm for r s0 = 0.1 cm to 0.06 cm for r s0 = 0.01 cm. 

In this work the fuel is chosen to be Hydrogen. This choice is made 
due to the relative simplicity of the Hydrogen reaction mechanism. The 
flammable mixture thus consist of dry air mixed with hydrogen, so the initial 
species are H 2 , 2 and N 2 . For all our simulations we take the equivalence 
ratio of 0.8. We have assumed that all the minor species in the air, such as 
argon and carbon-dioxide are negligible and that the initial gas is completely 
homogeneous. 

Radicals are not included in the initial mixture, even though the temper- 
ature profile already exists in the system at time zero. The radicals will start 
to form immediately after the simulation is started. 
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3. Model description 

In the following a new, novel, model for calculation of MIE is described. 
The model accounts for the chemistry only through the ignition delay time, 
and it only considers the central point in the temperature profile. Since a 
Gaussian temperature profile is assumed for all times, the spatial gradients 
are easily found when T m (t) and r s (t) are known. 

3.1. Obtaining an expression for the central temperature 

In this subsections it will be shown how the Gaussian temperature dis- 
tribution and the heat equation are used to find an expression of the central 
temperature as a function of time. In later subsections it will be illustrated 
how this expression is connected with the ignition delay time, and how the 
model determines a case of ignition. 

A large closed volume V is considered. Given that the initial temperature 
distribution is Gaussian it is assumed that the distribution stays Gaussian 
also during the short period until it can be determined if the mixture ignites 
or not. When the chemical reactions start influencing the temperature it 
is, however, clear that the profile will differ from Gaussian as the chemical 
heating initially will occur only in the center of the profile. The temperature 
distribution is then given by 

T (r, t) = (T m (t) - T ) e-(^) ? + T , (6) 

where To is the ambient temperature, T m (t) is the temperature in the middle 
of the distribution and r s (t) is the radius of the heat source, in this case the 
standard deviation. The initial values at t — are defined to be 

T m (0) = T max (7) 
r s (0)=r so . (8) 

The heat equation is required in order to include heat diffusion, and is 
given as 

§ = Q v>r, (9) 

where a is the thermal diffusivity. In 2D, it is convenient to use the polar 
coordinate system, and in 3D it is the spherical coordinate system which is 
most convenient. Since both 9 an are set to be symmetric, the derivatives 



5 



with respect to 9 and <p will be zero in our system. Hence, the Laplace 
operator is 

1 9 fr^f) , (10) 



where N = 1, 2 or 3 is the number of dimensions. 

By inserting the Gaussian temperature distribution, Eq. (J6|, into Eq. (J9|, 
and then evaluate the remaining equation at r = 0, it simplifies to 

~dT = -2«(T m )iV rg(t)2 . (11) 

Since a closed volume V is considered the specific volume n = V/M tota \, 
where M tota i is the total mass within the volume, must be conserved. This 
yields 

f 1 1 1Z f 

n — n = / dv = — / T — T dv = constant (12) 



v P Po P 

when n is the reference specific volume for T max = T and 1Z is the universal 
gas constant divided by the mean molar mass. Combining this with Eqs. (J6])- 
g gives 

(T m (t) - T ) • r s (tf = (T max - T ) • <. (13) 
From this equation r s (t) can be solved for, such that 

where again is the number of dimensions. The variable r s (t) can now 
be replaced in Eq. (11) to obtain a solvable first order differential equation, 
which is given as 

2+JV 

= -2a(T a )N ^- T "K . (15) 

1 s V- 1 max ^OJ 

The thermal diffusivity a{T) is here based on empirical thermodynamical 
data and fitted by the polynomial 

a (T m (t)) = (E + DT m (t) + CTl(t) + BT^(t) + Al£(t)) , (16) 

where the parameters are given in Table [1} 

To solve Eq. (15), which numerically is only zero dimensional, straight 
forward time stepping is used. The method applied in this work is the fourth 
order Runge-Kutta method. 
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Table 1: The parameters from the regression of a(T) 



Parameter 



Value 



A 
B 
C 
D 
E 



0.133608 
0.208675 
0.234699 
0.971533 
0.257533 



10 
10 
10 
10 
10 



13 



9 
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3.2. The ignition delay time 

The elements included in the model so far do not involve chemistry. To 
account for the chemistry involved in an ignition process the ignition delay 
time is used. 

To be able to predict ignition in a given mixture, one need to obtain data 
on the ignition delay time from that specific mixture. 

In order to measure the ignition delay time, a definition of when an ig- 
nition has taken place is required. There are several different definitions 
available in the literature, two definitions which are often used are 1) the 
time until the temperature has increased 400 K past the initial temperature 
and, 2) the time until the maximum temporal derivative of the temperature 
is achieved. 

Fig. [2] shows how the ignition delay time depends on temperature. In the 
figure there are two data sets, where both are obtained from zero dimensional 
simulations with the Pencil Code. The two are based on different methods of 
ignition determination. It can be seen that they differ when the temperature 
is sufficiently high. This is due to the fact that for high temperatures an 
increasing amount of the chemical energy is converted into radicals instead 
of thermal energy. This means that above a certain preheat temperature the 
mixture temperature will not increase with as much as 400 K, and this leads 
to the very rapid increase of r- m for definition 1). 

The model use r ig (T) together with the heat equation to predict if a mix- 
ture ignites or not, based on mixture composition and the initial parameters 
T max and r SQ . Lets define an ignition progress variable P such that P = 
at t = and P = 1 when the system has reached ignition. If P does not 
reach 1 this means that the mixture did not ignite, i.e. not enough energy 
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Figure 2: Ignition delay time as a function of temperature for two different definitions of 
the ignition delay time. The yellow line present the results found by defining ignition as 
the point where the temperature has reached 400 degrees above the initial temperature, 
while for the blue line ignition is defines as the time when the temperature gradient is at 
its maximum. The black line is a best fit to the blue line. 
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was supplied to the system. An equation, which fulfill these requirements, is 



F ^^» dt - (17) 



t=o 



The rationale behind this equation is that in order for ignition to occur 
a certain amount of radicals are required. These radicals are produced at 
all temperatures above a certain limit. The rate of radical production is 
inversionally proportional to the ignition delay time. 

Initially the parameters T max and r so are given. For each time step of 
lenght At, some heat diffuses away, and a new temperature, T m (t + At), 
is obtained from Eq. (TlSj) . Based on the new temperature a new ignition 
delay time, r ig (T m (t + At)), is obtained, which gives a new contribution to 
the ignition progress variable P. If, at any given time, P equals unity the 
current parameters and conditions have produced an ignition. As seen from 
Fig. [2j the lower the temperature gets, the higher the ignition delay time 



gets. This means that the term l/ri g (T m (t)) in Eq. (17) gets smaller, and 
each contribution towards P for each time step gets smaller. If heat is diffused 
away too quickly P will never reach 1, and the process will count as a non 
ignition case. 



4. Results 

In all of the following an hydrogen-air mixture with an equivalence ratio 
of 0.8 is considered. The ignition energy supplied to the mixture is 

Q= I c PP (T-T )dV' (18) 
Jv 

where V is a large volume containing the ignition source and T is given by 
Eq. ([5]). In order to find MIE for a given r s a series of simulations with 
gradually increasing T max are run. By definition MIE equals Q for the lowest 
temperature T max at which the mixture is ignited. In Fig.[3]the MIE is shown 
as a function of ignition source radius for 1, 2 and 3 dimensional setups. For 
the three dimensional case our DNS results with the Pencil Code (black line) 
is slightly above the results of Maas & Warnatz (1998) [3], this is however 
as expected since Maas & Warnatz (1998) jl] considered pure Oxygen as 
oxidizer, while in the current work air has been used. Furthermore, it is seen 
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Figure 3: MIE as a function of ignition source radius for different dimensionality's. Black 
line correspond to the DNS results from this work, blue line are the results for Maas & 
Warnatz (1998) [4] with a stoichiometric Hydrogen- Oxygen mixture while the red line is 
the results of Kim et al. (2004) 3J with a stoichiometric Hydrogen- Air mixture. 
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Figure 4: In this figure the results from the Pencil Code and the results from the simplified 
model are compared for all the dimensions. The parameters are To = 300 K, <fi = 0.8 and 
P = 1 atm. 



that the DNS results are very similar to the results of Kim et al. (2004) [3] , 
which was performed with a stoichiometric Hydrogen-Air mixture, except for 
the largest radius where the results of Kim et al. (2004) [3] bends upward. 
It is not known what cause the discrepancy for the largest radii. 

For the two dimensional results it is once again seen that the Pencil-Code 
results correspond very well with the results of Maas & Warnatz (1998) |1] 
except for the vertical shift due to the differences in oxidizer. 

For the one dimensional line there are no literature results with which it 
could be compared, but the trends nevertheless seems to be correct. 

Naively one could think that the MIE would scale as r^ where N is the 
dimensionality. Due to the stronger temperature gradients for smaller radii 
there will however be more thermal diffusion away from the center of the 
profile for smaller radii, and the MIE should therefore scale as where 
M < N. This is indeed also what is found in the simulations where M is 
0.83, 1.70 and 2.62, for ID, 2D and 3D, respectively. 

Lets define an ignition temperature T ign (r s0 ), which is a function of the 
initial width of the distribution r s0 , such that for T max < T ign (r s0 ) there is 
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Figure 5: Temperature evolution for DNS (black line) and model (red line) results. The 
time where the ignition progress variable P is unity is marked by the blue line. This 
example is taken from the one dimensional result with r s o = 0.03 cm. 



no ignition while for T max > T ign (r s0 ) the mixture ignites. In Fig. [5] T ign 
from the DNS simulations is shown as a function of r s o for the three different 
dimensionality's (black line). It is clearly seen that as the width of the 
initial profile is decreased the required temperature for ignition increases 
strongly. This is because the temperature diffusion rate scale as the second 
order gradient, which for a Gaussian distribution means that dT/dt ~ r s (t)~ 2 
(see Eq. (jTlj) ) . Furthermore it is seen that the higher the dimensionality the 
higher is the required initial temperature. This comes from the fact that 
high dimensionality yields more degrees of freedom for the thermal diffusion, 
and consequently the temperature in the middle of the distribution decreases 
faster. This is seen in Eq. (15) where dT/dt oc N. 

The red lines in Fig. [4] show the results obtained with the new model. It is 
seen that the model results are very similar to the results from the full DNS. 
Thus it is clear that for the setup used here the new model is a very good 
alternative to the full DNS, and there are no obvious reasons why this should 
not be the case also for other fuels, equivalence ratios, ambient temperatures 
etc. 
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In Fig. [5] the central temperature evolution is compared for the DNS 
results and the model results for the same specific case. The evolution is 
very similar for early times, but they are seen to deviate after ~ 1CT 4 s 
when the mixture starts to burn (as we see in the results of DNS) and thus 
produce heat. This is however not a problem since the model calculations 
will be finalized when P > 1, meaning that ignition has been achieved, 
which for this particular case happened after 9.5 x 10~ 5 s. If, however, the 
initial temperature is lower, such that there will be no ignition, the model 
calculation will be stopped, with the conclusion that no ignition could be 
achieved, when P has converged at a level below unity. 

5. Conclusion 

A new novel model for determining ignition has been described for the 
case of an initial Gaussian temperature distribution. The model, which re- 
quire a functional description of the ignition delay time as input, is compared 
against fully resolved DNS results and found to produce very similar results. 
This indicate that what determines if an initial ignition kernel will evolve 
into a successful ignition or not relies on two tings; 1) the amount of ther- 
mal diffusion away from the center of the distribution and 2) the integrated 
value of the inverse ignition delay. It is also expected that diffusion of radi- 
cals out of the center of the distribution will have an effect on ignition, but 
this is apparently only of equal or less importance compared to the thermal 
diffusion. 

Regarding the importance of the dimensionality of the heat source it is 
found that for higher dimensionality's there are more degrees of freedom for 
the thermal diffusion, and consequently the temperature will decrease more 
quickly in such a distribution such that the required maximum temperature 
for ignition is higher for the higher dimensions. 
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